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1 . INTRODUCTION 



The U.S. Navy has a continued interest in the reduction of the size and 
weight of propulsion components aboard both surface vessels and submarines. 
Some studies are currently underway in the Navy to determine the savings in 
both weight and volume which can occur if the condenser is designed to use 
corrugated or indented tubing to enhance heat transfer on both the shell and 
tube sides. While a comprehensive research program has been underway at the 
Naval Postgraduate School to study various enhancement techniques for 
horizontal condenser tubes (l,2] , an attractive alternative aboard submarines 
is a vertical condenser oriented aft of the steam turbine as shown in Figure 
l(a). With this arrangement, there would be a shorter machinery stack length 
compared to a horizontal condenser mounted underneath the turbine. This would 
allow a smaller diameter submarine to be built at a significant reduction in 
cost. The ultimate aim of this project is therefore to assess the reduction 
in condenser weight and size that may be feasible if the vertical tubes in 
this condenser are fluted on the outside (i.e., the steam side), Figure l(b). 

Many methods I 3 ] for enhancing condensation heat transfer have been 
proposed. Among them, Gregorlg [k] first recognized the importance of surface 
tension in film condensation on vertical fluted tubes. Thereafter, many 
experimental studies [5-ll] on vertical fluted surfaces were made to confirm 
his findings. Lustenader, et al., ( 5 ] employed a fluted vertical tube and 
obtained about four times larger heat transfer coefficients than those on a 
vertical smooth tube. Carnavos [6] carried out experiments on a doubly fluted 
vertical tube and found the augmentation ratio of condensation heat transfer 
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rate to be about six. Thomas [7,8] tested a vertical tube with longitudinal 
rectangular fins and wires. The augmentation ratio was about eight. Newson 
and Hodgson [ 9 ] manufactured various condenser tubes of highly enhanced heat 
transfer performance. Combs and his co-workers [10, 11 ] at Oak Ridge National 
Laboratory found the augmentation ratio to be around five for condensation of 
ammonia and refrigerants on vertical fluted tubes. 

Gregorig's theoretical model has also been improved upon but much still 
remains to be done. Edwards, et al., [ 12 ] proposed a condensation model on a 
heat transfer surface of triangular fins on the assumption of liquid film 
attaching to the tip of fin with a contact angle. The effect of a locally 
thin condensate film on the side of the fin (Figure 2) was not considered by 
them as well as by Fujii and Honda [ 13 ] . The latter, however, did consider 
two-dimensional conduction within the condensate film in the trough region 
(Figure 2) and within the fin wall. Mori, et al., [lU] considered the effect 
of a thin condensate film on the side of the fin but neglected the variation 
of its thickness in the vertical direction. They also considered only 
one-dimensional conduction within the fin and neglected any conduction within 
the condensate film in the trough region. Hirasawa, et al., [ 15 ] improved 
upon [ lU ] in that they included the variation of a thin condensate film in the 
vertical direction but neglected conduction within the fin and the film 
completely. Panchal and Bell [ 16 , 17 ] also neglected conduction within the 
fin and the film while analysing a sinusoidal fluted tube , but later found 
that two-dimensional conduction is important within the fin and the film for a 
triangular fin [ 18 ] . A recent empirical analysis by Barnes and Rohsenow [ 19 ] , 
based largely on (l 6 , 20 ] , reports an augmentation ratio of about fifteen for 
condensation of steam on a fluted surface while most earlier studies 
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considered condensation of a refrigerant and found much smaller augmentation 
ratios. With present day knowledge, it is difficult to ascertain whether this 
discrepancy is due to different fluids or due to a questionable analysis. 

All theoretical analyses discussed above break up the fluted surface into 
basically two parts. In the portion near the crest, gravity is neglected in 
comparison to the surface tension effect, while in the portion near the 
trough, the reverse is done. The patching between the two regions is carried 
out at a point that is selected quite arbitrarily at times [ 15-18] . This 
isolation of the two important effects, gravity and surface tension, is 
justified on the basis that the condensate film is thick in the trough region 
and thin over the crest. This is, however, not true over the initial portion 
of the tube length in the vertical direction nor during the initial portion of 
the tube just below a condensate drainage skirt (Figure 1(b)). These initial 
portions may be a few pitches in length. A recent analysis by Stack and 
Merkle [2l] does attempt to solve the complete equation with both 
gravitational and surface tension effects included over the entire flute but 
has three major drawbacks. First, when the condensate film in the trough 
region has thickened, the analysis treats it in the boundary- layer-sense , that 
is, it neglects all velocity gradients except those in the direction normal to 
the fluted surface. Second, the analysis is restricted to impracticably low 
values of amplitude to pitch ratio of the flute (0.02 and 0.0U) owing to the 
use of a Cartesian coordinate system rather than a surface-oriented coordinate 
system. Third, the analysis does not consider any heat transfer effects. 

Their analysis is only confined to finding the condensate-vapor interface. 
Fujii and Honda [ 13 ] also considered the entire thin film as one piece over an 
initial length of the tube (about 1/8 of the pitch) but neglected the 
important surface tension effect. 
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A major deficiency of all the above theoretical analyses except [21 ] lies 
in the way surface tension effect is determined. Contrary to the general 
belief, this effect does not depend upon the curvature of the condensate-vapor 
interface. Instead it depends upon the variation of this curvature along the 
interface. However, since the location of this interface is unknown a priori, 
many analyses [12, 13] simply determine it on the basis of the known flute 
shape owing to the argument that in the crest region where the surface tension 
effect is important, the film is very thin. Such analyses cannot be applied 
at all to triangular or rectangular fins since their curvature as well as the 
variation of curvature along the flute is zero. The same is, however, not 
true for the condensate-vapor interface. Even for sinusoidal flutes, there 
are large differences between the curvature and its variation along the curved 
surface for the given flute and the actual interface (see Fig. 10 and its 
discussion in Section 5). Moreover those analyses that do claim to determine 
the surface tension effect based on the actual shape of the condensate-vapor 
interface are also questionable in their claim as shown in Appendix A. 

Another problem with all earlier theoretical analyses except [ 13 ] is that 
conduction within the condensate film and the tube wall is either completely 
neglected or considered to be at most one-dimensional. Panchal and Bell [l8] 
point out clearly that two-dimensional conduction should be considered within 
both the condensate film and the tube wall. 

In an attempt to take care of these deficiencies, an analysis has been 
developed during this study that solves the complete equation with 
gravitational as well as surface tension effects included over the entire 
fluted surface using the surface-oriented coordinate system. In addition, for 
the first time, surface tension effects are determined as they should, from 
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the actual shape of the condensate-vapor interface rather than the shape of 



the flute. Two-dimensional conduction within the condensate film as well as 
the fluted tube is also included. A finite-difference solution of the highly 
non-linear partial differential equation for the condensate film thickness is 
coupled with a finite-element solution of the conduction problem. Details of 
the analysis follow in subsequent sections. 



7 



2. ANALYSIS 



Considering the condensate as a viscous, incompressible Newtonian fluid, 
and setting up a curvilinear orthogonal coordinate system (x^, x 2 , z) as shown 
in Fig. 3, the continuity equation yields [22) 



R + x„ 



3x, 



3u,. 

3x,. 



u. 



R + x. 



3w 

3z 



= 0 



( 2 . 1 ) 



where (u^ , u^, w) are the velocity components in the (x,, x^, z) directions 
respectively, and R(x^) is the radius of curvature of the fluted tube. The 
momentum equations are (cf. Ref. 22, p. 68) 



x^-direction : 



3u, 



u. 



R+x„ 1 3x 



3u, 



u 



2 3x, 



w 



3u^ 

nr 



U 1 U 2 

R+x^ 



+ — 



1 _R _3p_ 

p R+x 2 3x x 



= v 



3 2 u. 



(R+x 0 ) 2 3x 2 



3 2 u, 



3x,. 



3 2 u, 



3z 



R+x„ 



3x- 



u. 



(H+x 2 )‘ 



2R 



3u, 



(R+x 2 ) 2 3x l 



dR 






(R + x 2 ) 3dx l 2 



Rx 2 dR ^1 
( R+ x 2 )3 d X;L 3 Xl 



( 2 . 2 ) 



8 



X2-direction : 
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z-( vertical) direction : 
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where p and v are the density and kinematic viscosity of the condensate, p is 
the pressure in the condensate film, g is the acceleration due to gravity, and 
p v is the density of the vapor. 

In order to simplify these equations, we make the usual assumptions that 
inertia terms are negligible compared to other terms, and that 3/3x^ » 3/ 3x 
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or 3/9z, and u 2 << u, or w. This simplifies the momentum equations (2.2) - 
(2.4) to 





1 R 9p 
y R + x 2 9x 1 ’ 



(2.5) 




( 2 . 6 ) 





* 



(2.7) 



where y (=pv) is the dynamic viscosity of the condensate. Integration of 
(2.7) with boundary conditions 



yields 



9w 

9x„ 



w = 0 at x 2 = 0 , 



at x 2 = 6 



( 2 . 8 ) 
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Here S 3 is the (known) shear stress (in z-direction) on the condensate-vapor 
interface, and 5(x^,z) is the condensate film thickness. Similarly, 
integration of (2.5) with boundary conditions 



u j = 0 at x 2 = 0 



( 2 . 10 ) 



3u^ 

"3x1 
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y x 2 + 
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x 2 In 






R + x, 
R 



R + x. 



+ R In 



- x. 



( 2 . 11 ) 



Here is the (known) shear stress (in x ^-direction) on the condensate-vapor 
interface . 

It is advantageous to use the integral form of the continuity equation 
rather than the differential form (eqn. (2.1)). We therefore consider the 
control volume of Fig. h , and write the continuity equation as 





( 2 . 12 ) 



where kp is the thermal conductivity of the fluid, and hfg is the latent heat 
of vaporization. The right side of eqn. (2.12) represents the rate at which 
condensation of vapor takes place. The temperature gradient 3T/3 x 2 can be 
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approximated by (T s - T w )/6, where T s is the saturation temperature at which 

condensation takes place, and T w ( x^,z) is the wall temperature at x 0 = 0. 

This approximation basically assumes one-dimensional conduction in the 

condensate film but since 3/3x 2 >:> 3/3 x j or 3/3z, and we are interested only 

in (3T/3x_) at present, it is justifiable to make this approximation. 

2 6 

Substituting for w and u 1 from eqns. (2.9) and (2.1l) into eqn. (2.12), 
and integrating, we get 



3 1 35 * 
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2y 3x. 
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s -2 -2 3 S .3 

+ -1 1§- + i 1 £_ , ) ii. = _ 
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. (2.13) 



fg 



The pressure p can be related to the surface tension a and radius of curvature 
Rj_ of the condensate vapor interface by 



P = P s ± °/ R i > 



(2.1U) 



where p s is the saturation pressure of the vapor. The positive sign holds for 
0 < x^ < Xp/2, and the negative sign for Xp/2 < x ^ < Xp, where Xp is the 
length (curve DE in Fig. 3) along the fluted surface (in x ^-direction) from 
the crest to the trough (over half the pitch). We are basically interested in 
analyzing flute shapes that are symmetric about the crest and trough, and 
therefore take advantage of the symmetry by considering only half of the flute 
pitch. 
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Assuming o to be constant, the pressure gradient may be calculated from 



eqn. (2.14) 




It is thus the variation of curvature of the condensate-vapor interface that 
is of significance, not the curvature itself. In general, this can be quite 
different from d(l/R)/dx 1 , a fact neglected by many previous analyses. 
Relations for finding d(l/R^)/dx^ are given in Appendix A. 

Equation (2.13) is a partial differential equation for the condensate 
film thickness 6(x^,z). It requires the prior knowledge of shear stresses S, 
and S 3 that come from vapor dynamics and of T w (x^,z) that comes from heat 
conduction analysis for the fluted tube and condensate film. Eqn. (2.13) 
involves only first order derivatives w.r.t. z but fourth-order derivatives 
w.r.t. x^ owing to the presence of dp/dx^ (cf. eqn. (2.15) and (A-15) or 
(A-18)). For flute shapes that are symmetric about the crest and trough, the 
boundary conditions for the solution of eqn. (2.13) are 



86 8 3 6 




0 at x^ = 0, and at x^ 



x for all z; 
P 



(2.l6a) 



6(x 1 ,0) = 0. 

While the boundary condition at z = 0 is correct, it is not practical to start 
the integration of eqn. (2.13) from z = 0. Therefore, following Stack and 
Merkle l 21 ] , we replace the boundary condition at z = 0 by 

6(x 1 ,z (J ) = 6 Q , ( 2 . 16 b) 
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where <$ Q is the film thickness at the initial station z Q , and is found by the 
classical two-dimensional Nusselt solution 



Here T W g g — T w ( 0 , z g ) = T w ( 0 , 0 ) . 

For non-dimensionalization, we take the length L of the vertical fluted 
tube in the z-direction and length Xp and 6y as the characteristic lengths in 
the x 1 , and x 2 -directions respectively. Here 6y is related to L by the 
classical Nusselt relation (2.17) in exactly the same manner as 6 Q is related 
to z . Thus, we let 



With this non-dimensionalization, eqn. (2.13) and boundary conditions (2.16) . 

can be written as 




1/4 



(2.17) 



A = 6/6 



Y 



X = Xj/Xp , Z = z/L . 



(2.18) 




(2.19) 



3A 

3X 




0 at X = 0,1 for all Z 



(2.20a) 



1/4 




(2.20b) 



where 



14 



z 0 - V L ' 



D l = 



3a L 
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^ v p y 
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D = 6 y /x P ' 



Pc = P/xp , 




and 



l - for 0 < X < 1/2 
) + for 1/2 < X < 1 



Rn ~ R i/xp 

Equation (2.19) is a highly non-linear partial differential equation. It 
is solved numerically by the finite-difference technique as detailed below. 
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3 . FINITE-DIFFERENCE METHOD 



Equation (2.19) is parabolic in Z so that a forward marching scheme in 



the Z-direction can be used. Thus the choice of will affect only the 
region near z Q , and its effect will die out as z increases. This was indeed 
found to be true. Equation (2.19) as written is in conservative form. An 
equivalent but non-conservative form can be obtained by expanding the 



equation is still parabolic in Z but the mass is not identically conserved 
when the equation is integrated numerically [ 21 ] . 

We divide the interval 0 < X < 1 into n equal parts. Taking <5X and 5Z as 

the step sizes in X and Z-directions respectively, and using backward 

+ 

differencing in Z and mixed differencing in X-direction, we write the 
finite-difference form of eqn. (2.19) as 



derivatives in (2.19) and dividing by 3 A 2 . In that non-conservative form, the 




+ D. 
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+ E. 
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+ E. 
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(3.1) 
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where the subscripts i and k represent the location in X and Z directions 
respectively (cf. Fig. 5), and 3 is a number that is selected between 0 and 1; 
3=0 corresponds to backward differencing, 3 = 1 to forward differencing, and 
3 = 1/2 to central differencing in X. Since eqn. (2.20b) gives A T s for all X 
(i.e., for all i) at Z Q (say k = l) , we can march forward in the Z-direction. 
In other words, A f s at (k-l) for all i are known in eqn. (3*l). Thus (3.1) 
leads to a tridiagonal set of non-linear algebraic equations to be solved for 
(n+l) values of A f s at each location k of forward march in the Z-direction. 
This non-linear set is solved by linearization and successive iteration. Let 
us write eqn. (3.1) in the form 






A. 



\m 



1+1 ,k/ 



R. 

i 



(3.2) 



where ( ) m represents the values at the (current) mth iteration, and the 
coefficients A-j_ , Bj_ , , and depend upon the method of linearization. 

For linearization three methods were tried. In the first method, we 
s imply let 





m-1 




> 



J = 1,2,3, 



(3.3) 



where ( ) m_1 represents the known values at the (previous) (m-l)th iteration. 
The iteration is started by taking = A-j^-l for all i. With this 

linearization the coefficients in eqn. (3.2) are 
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and the boundary conditions in (2.20a) have been taken into account. This 
method works but has two drawbacks. Its convergence is slow, and it requires 
heavy under-relaxation for convergence. The latter implies that it is 
necessary to under-relax the values of A after every iteration according to 



U) m = ( A ) m_1 + X 



(A) m - (A) m_1 



(3.5)* 



where the relaxation factor X < 1. In fact, it was found that X decreases 
from about 0.5 at Z Q to very small values as Z increases. Typically X < 0.1 
for Z > 0.02. 

In the second method, we used Taylor series expansion in order to 
linearize terms containing A 2 and A 3 . Thus we let 







(3.6a) 



and 




m-1 m ? 

2 K,J (VJ - K,*) 



m-1 



(3.6b) 



With this linearization the coefficients in eqn. (3.2) are 
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* This equation is to be treated in the context of an equivalent FORTRAN 
statement . 
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\m-l 



(3.7b) 



(3.7c) 



(3.7d) 



(3.7e) 



(3.7f) 
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•3 P 3D (6Z) ? m-1 

R 0 " 4 0 ,k-l + S 3 4 0 ,k-l + 777 \m-l ~ E 1 6X * 1-2 4 M *1 ,k ) 






. ^ .m-lf , / 9E.. 9E \ 

( 4 0,k) [ 2 ( Vk ) + E l lx (l - 26) * e 3 + 52 ( IT + iz 2 )^. 



(3.7g) 



? ? 3D (6Z) 

R = A" 5 + E_ A , . + . 

n n,k-l 3 n,k-l )( / A 



- E lf (l - 26 >( 4 n-l,k)' 



4( A v ) r 
\ n ,k / 



. m-1 




* E i f (1 - 26) 



+• E^ + 5Z 




(3.7h) 



where G is again given by (3.^g) and the boundary conditions in (2.20a) have 
been accounted for . 

It is of interest to note from eqns. (3.7) that the Taylor series 
expansion was not used to linearize the right side of eqn. (3.1) as well as 
the terms involving Q on the left side of (3«l)» While such a linearization 
of the right side of (3.l) can be carried out without much difficulty, it was 
not found to be beneficial. The unnecessary complication was therefore 
avoided. The Taylor series linearization of terms involving Q in (3.1) was 
also tried but was found not to work at all. This method was found to be the 
best of the three. It allows a larger step size ( 5Z) in the Z-direction, and 
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a larger relaxation factor X than the first method. Overall, it is better 
than the first method by a factor of at least twenty, and was therefore 
adopted. The third method consisted in the use of Newton-Raphson method for 
solution of the non-linear equations. A general outline of this method is 
given in Appendix B. Though this method does not require any relaxation 
factor, it produced absurd results (A's oscilling with X after a few step 
sizes in the Z-direction ) , and was therefore abandoned. 
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4 . COMPUTATIONAL DETAILS 



Before the set of equations (3*2) resulting from eqn.(3«l) can be solved, 
we need to specify values of the dimensionless parameters E^, and D^. This 
requires the prior knowledge of shear stresses S,and S^, and the temperature 
T w (x^, z) of the condensate-wall interface. The computer code does have the 
provision for specification of shear stresses S, and but for the present 
results both these stresses were set to zero. This is because vapor dynamics 
that determines these stresses is beyond the scope of the present work. The 
computer code, however, does calculate T w (x 9 z) by considering two-dimensional 
conduction within the fluted tube as well as the condensate film. For this, 
the Laplace equation 



k 





(4.1 ) 



is solved subject to the boundary conditions 

3T 

= 0 at x = 0, P/2 where P = pitch of the flute, 

T = T s at the vapor-condensate interface, (4.2) 

and T = T c at the coolant-tube interface , 
where k is the thermal conductivity, and T c is the coolant temperature assumed 
to vary linearly from the coolant inlet to exit temperature. The last 
boundary condition in (4.2) needs modification due to a film resistance on the 
coolant side but such a specification either requires an ad-hoc value of the 
film resistance or a consideration of coolant dynamics that was also 
considered beyond the scope of the present analysis. Solution of eqns. ( 4 .1 ) 
and (4.2) was obtained by a finite element method in which the region 0ABC0 
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(Fig. 3) was divided into several linear triangular elements. Details for this 
can be found in any text on the finite element method. Some care is, however, 
required since the thermal conductivity of the fluid in region A3EDA is vastly 
different from that of the tube material in region DECOD (Fig. 3). The finite 
element method was preferred over the finite-difference method owing to the 
irregular shape of surface AB (Fig. 3). From this solution, we get the values 
of T w (x, ,z) on the surface DE (Fig.3). 

Since the solution of eqn. (3.1) requires T w (x^,z) which comes from the 
solution of (4.1) for which knowledge of 6-values is essential (for locating 
surface AB in Fig.3), the two solutions are coupled. Also, since the solution 
of (3.1) is found iteratively, it implies that eqn.(4.l) should be solved at 
every iteration. This places a rather prohibitive demand on the computer 
time. However, since (6Z) is small, of the order of 10“ 5 , changes in 6 are 
small at every iteration with the result that it is not really necessary to 
solve equation (4.1) at every iteration. Thus (4.1) was solved only once per 
step in the Z-direction. This saves considerable computer time since the 
finite element solution of (4.l) requires almost as much computer time as 
nearly one hundred iterations for solution of (3.1). Some sample runs were 
initially made to confirm that the error in 6-values due to this time saving 
feature was really negligible. 

4 .1 Condensate-Vapor Interface Profile 

As pointed out in Appendix A, one method for finding d(l/R^)/dXj is by 
use of equation (A. 8) which involves derivatives of f^(x), where f^(x) 
describes the condensate-vapor interface at any Z. Several methods were used 
to find these derivatives but none of them proved worthwhile since f^(x) is 
known at unequally- spaced values of x. The methods are listed below for the 
sake of completeness . 



24 



i) Cubic splines were fitted in a least squares manner through fj_ values 
to find f-j/ and f^". Then cubic splines were fitted through the previously 
computed f ^ values to find f^ 111 . A slight variation of this method was to 
fit cubic splines repeatedly through f^, fj/ and f^" values to find fj/ 11 . 

Both these methods gave different values for f^" and f-j/ f f pointing to their 
futility. 

ii) Parabolas were fitted through three consecutive points in succession 
on the condensate-vapor interface, and and dCl/R-^/dx^ were computed at the 
point central to each set of three points. This also led to inconsistent 
results . 

iii) An attempt to fit a truncated Fourier series (having 3-4 terms) 
through fi values was made. This was attempted since we were working with a 
sinusoidal flute. However, since the fit itself was poor, no attempt was made 
to calculate the derivatives. A better fit* (Fig. 6) was obtained by fitting 
( a l+ a^ cos 2ttx/P) through all (n+l) values of f^, where the coefficients a^ 
and a^ were found for a least squares fit. However, even this fit was 

considered inadequate to find d( l/R-j^ ) /dx ^ . A slight variation of this was to 

fit (a^+ a^ cos 2ttx/P) through the first n/2 values of f^ (for x < P/2) and a 
similar curve through the last n/2 values. This produced a very good fit on 
visual inspection (Fig. 7). However, using it to find d(l/R^)/dx 1 leads to 
wiggles in 6-values. Solution of equations (3.l) and (4.1) starting from Z Q 
leads to film shape rising suddenly in the trough region (Fig. 8) after a few 
(6Z) steps, thus violating the boundary conditions there. Starting from 
Z = 0.07, with results at Z = 0.07 found by better methods, the solution leads 
to such large wiggles that even negative values of 6 are found only a few 

steps downstream. In an effort to improve upon this situation, values of 6 

were computed to correspond to the fit in Fig. 7 before iterating. This, 
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however, slowed down the convergence so much that even two hundred iterations 
were insufficient for the same tolerance. In fact, it appeared as if 
convergence would never take place. Nevertheless, this exercise did point out 
the reason for the failure of the above technique, since the differences 
between the original 5-values and those conforming to the good ( ! ) fit in 
Fig. 7 were found to be large, as much as 20^. 

iv) Efforts to use lagrangian interpolation to compute f^ values at 
equi-distant x values, and then use techniques described in Section 4.2 to 
find f-j/, etc. also produced absurd results. Use of equation (A. 8) to 
compute d(l/Rj_)/dx 1 was therefore abandoned. 

It may be mentioned that while using eqn. (1.8) to find 
d(l/Rj_)/dx 1 , a futile attempt was also made to solve the partial differential 
equation (2.13) directly in terms of 5 3 rather than 5. 

4.2 Derivatives of Condensate Film Thickness 



Use of either eqn. (A. 15) or (A.l3) to find d(l/Rj_)/dx 1 requires the 
determination of first three derivatives of 5 w.r.t. x ^ . Fortunately 5 is 
known at equidistant values of x^. It is therefore easier to find these der- 
ivatives than those of fj_ w.r.t. x. Several methods for this determination 
were also tried with varying degrees of success . 

i) Derivatives of delta were found by fitting a cosine to the vapor- 
condensate interface in a manner analogous to that described in Sec. 4.1, but 
it resulted in absolutely useless results. 

ii) Another fruitless effort was to fit least squares cubic splines, 
repeatedly or otherwise, through 5-values so as to find 5*, 5" and 5 ? 1 1 . 
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iii) The first derivative of 5 w.r.t. x^ was found by solving the following 
set of linear equations [ 23 , p. 56 ] : 



t t t 

5. , + H<5. + 5. 
l-l 1 1+1 



_6 

5x 



-( 5 i+1 - | 5 i _ 1 j , i = 1 , 2 ,. . . ,(n-l) , 



(U.3) 



5 



o 



= 0 



1 



5=0 

n 



5x = x ( 6 X) , 

i. P 



and the second derivative, 5", by solving 



6 . , + 105. + 5.^ 
l-l 1 l+l 



12 



K) : 



( 5 i-l ‘ 26 i + 6 i+l)’ i= 1 . 2 ,...,(n-l). 



55 o 



It 



+ 6, 




(b.h) 



" " 12 

6 , + 55 = 

n -1 n 



( 5x l 



\2 



(Vi - \)- 
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Though these sets of linear equations are tridiagonal and can he solved 
easily, they do represent an increase in computer time. Besides the method 
yields rather poor results. 

iv) A method that gives better results than those obtained by above 
methods is based on passing least squares quartics through seven consecutive 
points at a time. Leaving the details to the reader (see also [2k] , p.l+92), 

t n f i f 

we present the equations for finding 6j_ and 6 ^ for i = o, 1 ,..., n. 

These equations make use of the boundary conditions (2.16a), and the fact that 

5 is known at equidistant x, values. 

▲ 

Relations for 6 f , accurate to o|fix j ^ , are 

«: = (22«._ 3 - 6 n us . 58 { ._ 1+ 58«. +1+ 6 U . t2 - 225 . +3 )/ 2526X3 , 

i = 3 j • • • j ( n-3 ) 5 

s' - 0 , s' = 0, 

o n 

6 1 = (- 586 o - 6 T 6 1 + 8o5 2 + 676 3 " 22 \) / 2 525x 1 , 

6 2 = (" 676 o “ 36 V 586 3 + 67 “ 226 5 )/ 2526x x , (U * 5) 

<-2 - ( 224 n -5 - 6TS n-k - 58S „-3 + 36 Vl + 6 TS n) / 2525lt l ' 
and 

Cl - ( 22S n-k - 6 ™n-l - 806 n -2 * 6lS n-l + *>«„)/ 252& 1 • 



-28- 



Relations for 6", accurate to o(5x 1 j 4 



are 



S 'i * (- 135 i-3* 6 TS i- 2 - 19 S i-l- 705 i - 195 i + l + 6TS i.2 - 135 i* 3 )' 



i = 3,4, 



s" = (-2 66 _+ 1346 - 386, - T06 )/ 1326x 2 , 

O \ J d 1 Of 1 



6^ = (-136^+ 6T5 3 - 326 2 - 35^ 19 6 q )/ 132 5x 2 , 



5 j = (-136 5 + 6T5 u - 196 3 - T06 2 - 326^ 67 6 q ) / 132 6x 2 , 



a" = ( -136 ,+ 6T6 n ,- 19 <5 T06 n „-326 ,+ 67 6 )/ 1326x 2 

\ n-y n-4 n-3 n-^ n-1 n/ 1 



n-1 



= (-135 ,+ 676 326 0 - 36 . - 196 )/ 1326x. 2 , 

\ n-4 n-3 n-2 n-1 n / 1 



and 



$" = ( -265 + 1346 - 386 “70 6 )/ 1326x 2 . 

n \ n-J n-d n-1 n/ 1 






132 6x 2 , 



• « • , ( n— 3 ) , 



(4.6) 
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Relations for 6 fff 



are 



f i f / 



_ + 6 . + 6 . — 6 . . — 6 . _ + 
3 i-2 l-l l+l i+2 



4+3 



/ 6 6x; 



i ~ 3 5 ^ 5 • ♦ • > ( n-3 ) ) 



Iff f ! t 
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6 6x^ 
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Use of these equations for finding 6' , 6" and 6' ' 1 



leads to mild 



fluctuations in 6-values after a few steps in the Z-direction. Their use had, 
therefore, to be abandoned. Less accurate relations involving an error of 
0 6 2 also yielded similar results. 

v) The best method for finding the derivatives of 6 for this problem 
turned out to be the use of central difference formulae 



<$! = ( 6 i+ i- 6 i_i)/ 2 6x i » i = l,2,...,(n-l) , 

(4.8) 



t T 




for the first derivative, and 

\ = ( 6 i+ i - 2 5 i + 6 i_i) / 6x i * i=l,2,...,(n-l) , 

5 o ■ 

S n = ( { „-l - S n)> 2 S \ • 

for the second derivative. The third derivative was found by applying (4.8) 
to 6" calculated from (4.9). This method gave the best results for 
5-distributions, and was therefore adopted. 
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5- RESULTS AND DISCUSSION 



The complete computer program was developed in two stages. In the first 
stage the solution of eqn. (3*l) alone was attempted. This program was tested 
for the vapor and tube properties for which results are available in [21 ] • As 
pointed out in "Introduction" , Stack and Merkle solve a much simplified form 
of (3.1) since besides other simplifications, they use Cartesian coordinates 
rather than curvilinear coordinates. In fact, we had modified our earlier 
program to solve Stack and Merkle' s equation (l4) for the condensate film 
thickness, and obtained identical results. The finite element program for the 
solution of eqn. (4.1) was then developed, and the two were coupled. 

In the following we present results for one vapor-tube combination for 
which experimental data is also available 1 10 ,11 1 for comparison. Since our 
analysis assumes the radius of curvature of the fluted tube to be large 
compared to the film thickness, it is not realistic to consider sinusoidal 
flutes with amplitude-to-pitch ratio much above 0.2. With this in mind, we 
computed results for tube 'F' of Ref. [ 10 , 11 ] - This tube has seven skirts and 
an I.D. = 22.9 mm. Details of the actual shape of the flute, and its 
dimensions are not available in [ 10 ,11 ] , and were therefore approximated from 
a blow-up of the tube 'F* photograph in [ 11 ] . The relevant fluid (R-113) and 
tube (aluminum) properties are 



P = 1498.343 kg/m 3 , 

v = 3. 2067x10- 7 m 2 /s 
h fg = 145225-56 J/kg , 
k f = O.O6951 W/m-K 



^woo 



= 318.5 K 



p v = 8.58628 kg/m 3 , 

c = 0.01432 N/m , 

T s = 325.5 K 
k w = 205 W/m-K 

T c = 318.5 K (in), 318.8 K (out) , 
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L = 0.1H2875 ra | = x ) , h t = 0.88077 mm , 

P = 1.6lUl mm , h,Q — 0.155^3 mm 5 

where kf and k w are the thermal conductivities of the condensate and fluted 
tube respectively, and h^ and h,-, are associated with the flute shape (Fig. 3) 
taken to be 



f(x) = h^ + h Q cos . (5.1) 

Thus hQ is half the amplitude of the flute. The fluid properties were 
calculated at 322 K (mean of T s and T woo ) from known correlations ( 25 ] • These 
values lead to 



6y = 0.08051 mm , x 0 = O. 876 H 87 ram 

After some numerical experimentation involving different step sizes, 
etc., the solution of eqn. (3.1) was started from Z Q = 5x10“ 6 with g = 1/2 
(corresponding to central differencing in X), 6X = 0.05, an initial <5Z = 
5x10“^, and an initial A = 0.5 (see eqn. (3.5) )• Further iteration for the 
solution of (3.1) was terminated when 



1 




< e for all i at every k , 



(5.2) 



where e was taken to be 10 - ^. As we marched downstream in the Z-direction, 5Z 
was increased and the relaxation factor A had to be decreased according to 
Table 1, which also gives an idea of the number of iterations before (5*2) was 
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satisfied. Clearly, it takes longer to get the converged solution as we 



proceed downstream in the Z-direction. Reasons for this will be apparent 
soon . 



Table 


1 


Some Parameters 


for Solution 


of 


(3.1) 


5Z 


X 




z 


ff 


of iterations 


5xlO" 6 


0.5 


up to 


0.002 




~ 15 


8xl0- 6 


0.2 


up to 


0.01 




~ 30 


1.25xl0- 5 


0.1 


up to 


0.0225 




- 50 


2xl0“ 5 


0.06 


up to 


0.086 




~ 90 



Three divisions of the region OABCO (Fig. 3) into finite elements for the 
solution of eqn. (U.l) were tried. The one selected on the basis of adequate 
computational accuracy and relatively economical calculation had 15^ 
triangular elements with 108 nodes. There were 11 nodes each on the faces AB, 
BC, OC and DE (Fig. 3), and 15 nodes on the face AO. 

Fig. 9 shows, to true scale, the flute shape (thick curve), and the 
condensate film shape (dashed curves) at four values of Z = 0.001, 0.01, 0.03, 
and 0.086. As expected the film thickens quite rapidly in the trough region 
while it remains thin over the crest. Unfortunately we are unable to present 
any results for Z much greater than 0.086 since convergence of our solution 
for the film thickness (with the same tolerance e) becomes very slow. The 
principal reason for this slow convergence is the thickening of the film in 
the trough. Moreover, it is really futile to try to speed-up the convergence 






since our formulation for the condensate film thickness is on shaky ground as 
the film thickens. The non-linear partial differential equation (2.13) comes 
from a houndary-layer type analysis, i.e., we neglect all velocity gradients 
except those in the direction normal to the fluted surface. While this is 
perfectly reasonable over the initial portion of the tube height when the 
condensate film is thin, it is increasingly questionable as the film thickens 
in the trough region. It appears that for the fluid-tube combination 
considered here, it is improper to analyze the film in the trough region in 
the boundary-layer-sense for Z > O.O 85 . It is also clear from Fig. 10, that 
shows the absolute value of d( l/R-j_) /dx^ at various Z values, that this 
derivative is negligible in the trough region as compared to its value in the 
crest region. It is therefore appropriate to neglect the surface tension 
effects in the trough region once the condensate film has thickened. Fig. 10 
also plots the absolute value of the derivative of flute curvature with X]_, 
and even at the low value of Z = 0.002, it is very different from the 
derivative of condensate-vapor-interface-curvature with X]_. 

While these considerations will be undertaken in an extension of this 
project, it is encouraging to note from Fig. 11 that a linear extrapolation, 
on a log-log basis, of our results to date shows a fairly good comparison with 
experimental data. Fig. 11 shows the heat load per half flute (in W) as a 
function of Z on a log-log scale. The solid portion of the straight line is 
based on computed results while the dashed portion is the extrapolation. The 
lone circled point is based on the experimental data ( 11 ] . We also checked 
that heat transfer rates across the faces AB and 0C (Fig. 3) matched within 
1 %. Theoretically they should be identical since faces 0A and BC are 
insulated. 
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6. SUGGESTIONS FOR FURTHER WORK 



Besides the deficiencies mentioned above, our present analysis also 
neglects the transverse curvature term, i.e., it assumes the film thickness to 
be small compared to the radius of curvature of the flute. The analysis is 
therefore -restricted, in its present form, to flutes having a small amplitude 
to pitch ratio. This is, however, a minor drawback since it seems 
(analytically at least) that we can overcome this deficiency relatively 
easily. Its implementation in the computer code is yet to be done. 

Another assumption, common to all analyses to date, and far more serious 
in terms of correct modelling of the practical applications, is the complete 
neglect of vapor shear on the interface. Our code does have the provision for 
studying its effect but since we have not analyzed the vapor dynamics yet, we 
simply take the vapor shear to be zero at present. One should also consider 
the coolant dynamics for a complete solution but all this will undoubtedly be 



very demanding. 
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APPENDIX A 



CURVATURE OF CONDENSATE- VAPOR INTERFACE 
AND ITS VARIATION ALONG FLUTE SHAPE 

Equation (2.13) for the determination of condensate film thickness 

6(x_^,z) involves the value of pressure gradient dp/dx which, in turn, 

depends upon d(l/R )/dx , as is clear from (2.15)* We present below some 

i 1 

relations for calculation of this quantity. 

The basic terminology is present in Fig. 3, wherein, f(x) denotes the 
flute shape and f^(x) the condensate- vapor interface. The position vector to 
a point on the flute is r w , and that .to a point on the condensate-vapor 
interface is r^_ , such that 

f r y — n w 5 9 ( A • 1 ) 

•* 

where n w is the unit normal to the flute as shown in Fig. 3. We also let Xj_ 
be the curvilinear coordinate along the condensate-vapor interface. 

Then by definition 



1 . 32? i 

*i " ax. z 


(A. 2) 


1 








Now r^ = (x, fp , 

1/2 

x » p 

and ^ (1 + fi ) d£ , 


(A. 3) 


(A.U) 



where prime denotes differentiation with x. From equation (A.U), we get 



dx-: 



dx 



(i + fl 2 ) 



L /2 



(A. 5) 
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From (A. 3) then 



dr a 



dr dx 



dx.: dx dXi 



, 2 - 1 /2 

= (1 + fi ) (l,fi ) 



Thus 






2 2 
dx i (l + fi 2 ) 



(-fi , 1) 



(a. 6) 



From equation (A.2) then 



R i 



(1 + fi' 2 ) 



3 / 2 



And 



(A. 7) 



d (i-) - 4- (i-) dx 



dx]_ R^ dx dX]_ 



sgn( f j_ ) 



(l+fi ^ 



T/2 



. "2 ”1 



fi - 3 



f . f . 



(l+fi' 2 ) 



(1+f 2 ) 



— *71 ’ 

io. 1 ^ 



(A.8) 



dx 

]_ '2 */2 

since ,j x = (l+f ) in analogy with (A.5)» 



(A. 9) 
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While equation (A. 8) seems to he in a convenient form for calculation of 

d(l/R-;)/dx , it is not practical since f,- is knovn only at a discrete set of 
1 1 1 

f TT t? T 

unequally-spaced x-values . Thus , finding fj_ , f ^ and f ^ is troublesome , to 
say the least. Several attempts to use equation (A. 8) had to be abandoned for 
this reason (see Sec. H.l for details). 

Another relation for the curvature of the condensate-vapor interface 

-*• 

comes from use of equation (A.l), and the fact that the unit normal vector n w 
is given by 



( df dx n 
nw " dx L ’ dx£ ) 



_l/2 

= (1+f ’ 2 ) (-f ' , 1) . 



(A. 10) 



Thus, from (A.l) 



r i = x - 



6f 



d + f' 2 > 



T/2 • f + 

TO l *- 



’2 1/2 
(1+f 2 ) 



, since r w = (x,f) 



dr 



Therefore , 



dx 



1 / , ,» . 

= ( 1 - 6 f - 



6f 



(1+f' 2 ) 



V 2 



f'+ 6* - 



T TT 

5f f 



(1+f d ) 



Tfl 



(A. 11) 



where 5 = 86/3x_^ but f = df/dx, etc.. The reason for keeping 36/3x_^ in place 
of 36/3x is that solution of equation (2.13) leads to values of 6 at equidis- 
tant x_^-values but at non-equidistant x-values. It is therefore convenient 
to find numerically derivatives of 6 with respect to x . 
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Carrying on further, we can find d 2 r.j_/dx 2 from (A.ll), and then use the 
relation 



1_ 

R,- = 



Vl 2 l?iT - (?i:?i'’) 2 



l /2 



Ta 



(A. 12) 



where 




drj_/dx, to get after some lengthy algebra the expression: 




Y^S + y ^S 6 + + Y^ (l-66 + 26 ^ )| 



[yg(l+6 2 ) + y^S 2 + Y g<$] 



3 / 2 



(A. 13) 



where 




Y p = (—V) 
1+f 2 



Y 3 = 



-2f 



"2 



(1+f 2 ) 



3 /2 






-o " 1/2 

= (1+f 2 ) (f 



3f f 



1+f 



Y 



5 



(1+f' 2 ) 



3 / 2 





y 8 



2f 



(1+f 2 ) 



T /2 ‘ 
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It is interesting to note that many theoretical analyses that claim to 
determine the surface tension effect based upon the actual shape of the 
condensate-vapor interface invariably use 




d+6' 2 ) 



V: 



(A.lU) 



for the curvature of the interface. This follows from the correct relation 
(A. 13) if only the first terms in both the numerator and the denominator are 
retained. However, even for small 5, there is hardly any justification for 
dropping all other terms in (A.13)« The difference between the actual 
d(l/Rj_)/dx and that obtained from (A.lU) gets further amplified due to 
differentiation . 

For convenience let us write 




N.sgn(N) 
/ 2 



where the expressions for N and D are clear from equation (A. 13). It is then 
easy to get 



d ; 1 i _ sgn(N) 

dx T Tl 



5"'(y 5 -5y 1 ) + <$"{Y 5 ' + + 3Y-L5' } + \ 



+ 6 2 (y 1+ + 2y 1 ') + 6' (y 3 + S(2y 2 + Y^ ' ) ) + S 2 Y 2 ' + 6 Y 3 ' 



- 2rTd" {S'(2YgS" + 2 y 7 <$ + Yq) + (l + 6 2 )Yg' + S 2 Y j' + SYg'}, (A.15) 

I 
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Where 



and 



dx 



dYi , - l h 

= f (l+f 2 ) 



. dY 2 3f" 2 



'2 dx. 



,» i 2f f 



T _T»2 



1 (l + f’2) 7/2L 1+f2j 



, d ^3 

r 3 “ dx. 



2f 



, 2 

(l+f 2 ) 



1 Ito 

3f f 2 



l+f 2 



- 2f 



t dy U fiv f" . up i in. 9f'2f"3 

y U = dF" = ~ ± ~ T ~ ~ (3f 2 + 7f f ) + 



11 « 2 ( l+f ' 2 ) 2 



(l+f' 2 ) 3 



dYc 

! S t tT 

Y 5 = df =3ff 



dY, 



Y g «^-2f'f"(l+f' 2 ) 



io-V2 



dY, 



1 - _2f 



'7 dx. 



(l+f c ) 



M t »T2 

, 2f f i 

5/2 (f '2 

'2x 7 l+f 



dYc 



t _ii 2 



f f 



8 “ dx 



it i , 

- f ) 



1 l+f l+f 
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Iquation (A. 15) is quite complicated but as mentioned in the discussion on 

3 3 

qn. (2.13), it contains 8 5/8x]_ . It is therefore essential to determine the 
ondensate film thickness as accurately as possible. Otherwise, errors in 
inding the third derivative of 5 w.y*t. X]_ numerically will be excessive, 
ue to these numerical difficulties, it was found that use of eqn. (A. 15) to 
ind d(l/Ri)/dX]_ resulted in mild oscillations in 5-values after a few steps 
n the Z-direction. Several methods (at least five) were used to find 
-derivatives w.r.t. X]_ but all resulted in oscillatory 5-values; the least 
scillations resulted from the use of central-difference formulae. 

In order to overcome this difficulty, another relation for 1/R^ and its 
erivative w.r.t. X]_ was developed. This relation, based on eqns. (A. 2) and 
A.l ) , is 




3 2 (? w + 



i^S) 



3x,- 



3 2 ? w 3 2 n w 8n w 

+ 6 + 2 

3x^ 2 3x^ 2 3Xi 



36 



“v 



3 2 6 



5 *i 2 



(A. 16) 



here 3/3x^ has heen approximated by 3/3xj_. Noting that r w = (x,f), and riy. 
3 given by eqn. (A.IO), we get from (A.l6), after some algebra, that 



" "I 1/ 2 

= ( 0 ^+ a 2 6 + a 3 6' + c^fi") 2 + ( ^ + ^6 + & 3 6* + g^d") 2 , (A. 17) 
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where 6' = 36/3x , f' = df/dx, etc., and 

1 



a = -8 f 

1 1 

“2 = 7T7I 

(1+f 2 ) 



a 3 = -2B 1 



a U = f' , 



M . ’a "2 

8i = f (1+f 2) 



S>2 ~ a 2^ “ Bi^f , 



UfY’2 ... 



— - f ) , 



i + f 



B 3 = 2ct]_ , 



and 



84 = (1+f 2 ) 



~ l /2 



From equation (A. 17) it follows that 
^~(p7— ) = + a 2 5 + a^6 + 6 ) ( 8 2 + ^ ' 5 + — 3 8^ 6 + 

+ (8^ + + 8^6 + 8^ 6 )(-a 2 + 8j 6 + 382^ + 3 a ^5 + 8^5 )} 



»T f 

V ) 

(A.l8) 



where 



a,. 



da 2 

dxn 



_1 f'Uf'g + 13f ’ f " ’ ) _ fiv _ 28 f'2f"3 



( 1+f ' 2 V 



1 + f 2 



and 



d0 2 

0 = ^ = a f + f (20 a - a 0 ) - 0 0 f 

2 ^ — 2 U 2 11 1 h 

1 



Equation (A.l8) also contains 3 6/3x , but unlike equation (A. 15), it yields 

stable (non-oscillatory ) values of 6 upon integration of eqn.(2.13). This 
relation was therefore used in the final computer code. 
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APPENDIX B 



Nonlinear Equations Solution by Newton-Raphson Method 



Consider the set of equations 



[K] {x} = {R} 



n 



or E kj_j(x , x^, x n )xj = R^x^ , x^, x n ) , i = 1, 2, n. 



j=l 



1 2 



where k,- t and R,- are non-linear functions of x , x , . . , x,_ . 

- L J 1 2 11 

n 

Let F i ( x , x > . . , x n ) = E k.! 1 ( x , x , . . , x ) x i - R,- ( x , x , . . , x n ) , 

1 1 2 n j =1 1J 12 n J 1 1 2 n 

l — 1) 2 , n « 



th 

So Fj_ is the residue in the i eqn. For solution Fj_ = 0, i = 1, 2, n. 
Expanding Fj_ by Taylor series , we write 



F,-(x + Ax , x + Ax , x n + Ax n ) 

i 1 q 2 2 n n 



Fj( x , x , 
1 1 2 



, x n ) + E 9F i (x , x , 

0 = 1^1 1 2 



) Ax< 



i = 1, 2, 



• II • 



Neglecting all higher-order terms and setting leftside = 0 for a solution 
we get 
n 



£ 8F i (x , x , 
j=l 9xj 1 2 



x n ) Ax * = -F i ( x , x . . . , 
n j 1 2 



x n ) 



i = 1 , 2 , . . , n. 
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Thus Newton-Raphson iteration algorithm can be written as 



, , m m+1 m 

[J] {Ax} = -{F} 



m+1 m m+1 

{x} = {x} + {Ax} 



, m 



where m denotes m™ iteration & [j] is the system Jacobian at the m™ itera- 



tion. 



m 



The elements of [j] are 



m 

J ij 



3F 



3x 



- (x? , x 2 m , .. , x£) . 



m 



If it is too time consuming to compute [j] at every iteration, find [ J ] once 
and use it until convergence, which will be slow, however. 
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SINGLE 

CONDENSER 



(a) Condenser arrangement 




(b) Cross section of vertical condenser. 
Fig. 1. Vertical Submarine Condenser 
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CONDENSING SIDE 
(fluted outside surface) 
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Fig. 2. Fluted tube principle of operation (condensation mode) 
surface tension forces acting to push condensate from crests into troughs. 



y 




Fig. 3 Co-ordinate System. 
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Fig. 4 Control Volume in the Condensate Film. 
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FIG. 5 Finite Difference Grid. 



53 




(ujuj) X 
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Fig. 6 One Cosine Fit (---) Over the Entire Condensate- 
Vapor Interface ( — ) at Z — 0.08. 
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Two Cosine Fits ( — ) Over the Condensate -Vapor 
Interface ( — ) at 0.08. 




Flute Shape ( — ) and an Erroneous Film Shape 
at Some Z . 



y (mm) 



1.2 




x (mm) 



Fig. 9 Flute Shape (— ) and Correct Film Shapes 
at Some Z Locations . 
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Fig. 10 Comparison of Dimensionless Derivative of the Curvature of 
Flute ( — ) and Condensate -Vapor Interface at Various Z Values. 




Fig. II. Heat Load(inW) as a Function of Z Computed j , 

Extrapolated , Experimental [ll_ 
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